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Abstract 


We  have  developed  an  intermolecular  potential  that  describes  the  structure  of  the  a-form  of 
the  hexahydro-l,3»5-trinitro,l,3,5-s-triazine  (RDX)  crystal  The  potential  is  composed  of 
pairwise  atom-atom  (6-exp)  Buckingham  interactions  and  charge-charge  interactions.  The 
parameters  of  the  Buckingham  repulsion-dispersion  terms  have  been  determined  through  a 
combination  of  nonlinear  least-squares  fitting  to  observed  crystal  structures  and  lattice  energies 
and  trial-and-error  adjustment.  Crystal-packing  calculations  were  performed  to  determine  the 
equilibrium  crystallographic  structure  and  lattice  energy  of  the  model  There  are  no  significant 
differences  in  the  geometrical  structures  and  crystal  energies  resulting  firom  minimization  of  the 
lattice  energy  with  and  without  symmetry  constraints.  Further  testing  of  the  intermolecular 
potential  has  been  done  by  performing  symmetry-constrained  isothermal-isobaric  Monte  Carlo 
simulations.  The  properties  of  the  crystal  (lattice  dimensions,  molecular  orientation,  and  lattice 
energy)  determined  from  Monte  Carlo  simulations  at  temperatures  over  the  range  4.2-300  K 
indicate  good  agreement  with  experimental  data.  The  intermolecular  potential  was  also  subjected 
to  isothermal-isobaric  molecular  dynamics  calculations  at  ambient  pressure  for  ten^eratures 
ranging  from  4.2  to  325  K.  Crystal  structures  at  300  K  are  in  outstanding  agreement  with 
experiment  (within  2%  of  lattice  dimensions,  and  almost  no  rotational  and  translational  disorder 
of  the  molecules  in  the  unit  ceU).  The  space-group  symmetry  was  maintained  throughout  the 
simulations.  Thermal  expansion  coefficients  were  determined  for  the  model  and  are  in  reasonable 
accord  with  experiment. 
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1.  INTRODUCTION 


The  pursuit  of  determining  the  microscopic  details  of  the  chemical  and  physical  processes 
that  occur  in  condensed  phase  energetic  materials  has  led  to  the  development  of  a  potential 
energy  function  to  describe  the  a-form  of  the  hexahydro-l,3»5-trmitro-l,3,5-j-triazine  (RDX) 
crystal.  Crystalline  RDX  exists  in  two  phases  [1]:  the  ambient  phase  (a-solid),  for  which  the 
structure  has  been  characterized  by  neutron  diffraction  measurements  [2],  and  an  unstable  phase 
(P-solid),  the  crystal  structure  of  which  has  not  been  determined.  RDX  is  one  of  the  most  widely 
used  explosives,  and  numerous  data  on  the  physical  and  chemical  properties  exist  for  this 
material  [3].  However,  direct  measurements  of  molecular-level  details  of  the  response  of  RDX 
to  external  stimuli  (such  as  heating  or  shock)  are  not  available.  Simulations  using  models  such 
as  that  presented  here  provide  invaluable  insight  into  the  dynamic  responses  of  the  energetic 
material.  This  information  can  be  used  to  determine  how  to  manipulate  or  control  the  behavior 
of  the  material.  Toward  this  end,  we  present  an  atomistic  model  for  the  intermolecular 
interactions  in  the  a-RDX  crystal  [2].  We  will  describe  the  procedure  used  to  develop  and 
parameterize  a  simple  functional  form  of  the  intermolecular  potential  energy  in  the  first  part  of 
this  paper,  and  then  describe  a  series  of  tests  to  which  this  model  was  subjected  to  determine 
whether  it  accurately  represents  the  a-RDX  crystal. 

Theoretical  studies  of  the  past  two  decades  have  demonstrated  that  simple  atom-atom  pair 
potentials  are  sufficient  for  predicting  molecular  crystal  structures,  but  the  proper  development  of 
accurate  models  is  challenging  [4].  The  main  difficulty  lies  in  attaining  a  proper 
parameterization  of  the  potential  function  such  that  the  model  reproduces  properties  of  ihe  crystal 
such  as  the  lattice  energy,  crystal  structure,  elastic  properties,  and  phonon  frequencies.  Often,  the 
resulting  potentials  can  be  used  to  predict  properties  not  included  in  the  original  fitting  of  the 
function.  Also,  when  large  numbers  of  crystals  with  similar  functional  groups  are  used  in  the 
fitting,  the  potential  parameters  are  transferable  [4-6].  Reproduction  of  properties  used  in  fitting, 
prediction  of  properties,  and  transferability  of  potential  parameters  are  three  metrics  used  to 
assess  the  quality  of  a  model  potential. 
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Numerous  studies  of  organic  molecular  crystals  [4]  have  demonstrated  that  in  many  cases,  the 
intermolecular  interactions  can  be  described  by  using  simple  isotropic  potentials,  such  as  the 
(6-exp)  Buckingham  potential 


V  =  Aexp(-Br)-C/r^ 


(1) 


or  the  (6-m)  Lennard-Jones  potentials 


V  =  D/r’“-C/r^ 


(2) 


where  m  =  8-14. 

The  repulsion  and  dispersion  parameters  A,  B,  C,  and  D  are  assumed  to  depend  only  on  atom 
type  and  thus  are  transferable.  Despite  their  analytical  simplicity,  these  potentials  have 
accurately  described  a  large  class  of  crystals,  particularly  nonpolar  crystals.  A  number  of  sets  of 
such  empirical  intermolecular  potentials  for  organic  crystals  containing  C,  H,  N,  O,  Cl,  and  S  are 
currently  available  [4, 6]. 

A  different  approach  for  describing  molecular  crystals  is  to  consider  explicitly  the 
electrostatic  interactions  between  the  atoms  of  the  molecules.  In  the  simplest  case,  the  exp  6-1 
potentials 


V  =  Aexp(-Br)-C/r^  +  E/r,  (3) 

or  m-6-1  potentials  offer  a  significant  increase  in  the  flexibility  and  transferability  of  the 
nonbonded  potential  parameters.  This  form  includes  the  Coulombic  interaction  between  the 
charges  associated  with  the  various  atoms  in  a  molecular  crystal.  As  an  example,  Williams  and 
Starr  [7]  have  shown  that  the  (exp-6- 1)  potential  gives  excellent  results  for  calculations  of  lattice 
vibrational  for  aromatic  hydrocarbons,  and  many  of  the  available  force  fields  such  as  Amber  [8], 
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ECEPP  [9],  or  Dreiding  [10]  use  this  type  of  potential  term  for  simulations  of  organic,  biological, 
and  main-group  inorganic  crystals. 

Despite  the  success  of  this  kind  of  representation,  many  crystals  contain  substantial 
anisotropies  in  the  interactions  between  the  molecules  in  the  crystal.  Such  crystals  cannot  be 
adequately  described  using  the  simple  electrostatic  multipoles.  Significant  progress  has  been 
made  by  the  use  of  the  distributed-multipole  technique,  which  has  been  used  to  successfully 
describe  hydrogen-bonded  structures  [11]  and  tt-k  interactions  [12],  A  recently  proposed 
algorithm  for  relaxation  of  molecular  crystal  stractures  using  a  distributed  multipole  electrostatic 
model  provides  a  promising  approach  for  studying  packing  of  polar  crystals  [13], 

The  main  objectives  of  the  work  presented  here  were  to  develop  an  interaction  potential  of 
the  a-RDX  crystal  that  both  accurately  reproduces  experimental  information  and  is  simple. 

These  objectives  must  be  accomplished  since  the  aim  is  to  use  this  model  in  predicting 
phenomena  occurring  in  condensed-phase  RDX,  which  might  be  difficult  to  measure 
experimentally.  If  the  model  is  too  complex,  CPU  requirements  to  calculate  bulk  properties  of 
the  crystal  would  be  prohibitive.  If  the  potential  energy  function  is  too  simple  and  cannot 
reproduce  measured  data,  then  the  model  is  not  representative  of  RDX  and  thus  is  not  useful. 

We  have  found  that  the  RDX  crystal  can  be  reasonably  described  using  the  exp-6-1  potential 
(equation  [3]).  Several  of  the  parameters  (such  as  the  Coulombic  terms)  were  assigned  before 
the  fitting  of  equation  (3)  to  experimental  observables  was  attempted.  The  Coulombic  terms 
were  determined  through  fitting  of  partial  charges  centered  on  each  atom  of  the  RDX  molecule  to 
a  quantum  mechanically  derived  electrostatic  potential.  We  also  used  previously  published 
parameters  for  H-H  and  C-C  interactions,  and  assumed  traditional  combination  rules  to  obtain 
heteroatom  parameters  from  the  homoatom  parameters.  The  remaining  parameters  were  selected 
such  that  symmetry-constrained  molecular  packing  calculations  reproduced  both  the 
crystallographic  structure  and  the  lattice  energy  of  the  crystal.  Throughout  all  calculations 
described  hereafter,  the  molecules  are  assumed  to  be  rigid,  and  the  structure  is  described  by  die 
center-of-mass  positions  and  orientational  parameters  for  each  molecule  in  the  unit  cell.  The 
intermolecular  potential  was  subjected  to  additional  tests:  molecular  packing  calculations 
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performed  without  symmetry  constraints,  symmetry-restricted  isothermal-isobaric  Monte  Carlo 
calculations  (NPT-MC),  and  isothermal-isobaric  molecular  dynamics  calculations  (NPT-MD). 

hi  section  2,  we  give  details  of  the  analytical  form  of  the  intermolecular  potential  used  in  the 
calculations,  the  accelerated  convergence  technique  employed  for  evaluation  of  the  lattice  sums, 
and  parameterization  of  the  potential  function.  Section  3  is  devoted  to  analysis  of  the  theoretical 
methods  used  for  minimization  of  the  lattice  energy  and  the  details  of  the  NPT-MC  and  NPT- 
MD  calculations.  Results  of  crystal-packing,  NPT-MC,  and  NPT-MD  calculations  are  given  in 
section  5.  The  conclusions  drawn  from  the  study  are  given  in  section  6. 

2.  INTERMOLECULAR  POTENTIAL 

The  central  problem  in  classical  simulations  of  molecular  crystals  is  the  construction  of  an 
analytical  potential  energy  function  that  accurately  represents  the  intermolecular  interactions  and 
thus  predicts  the  stractural  and  thermochemical  parameters  of  the  crystal.  In  this  work,  we  adopt 
some  general  principles  for  atom-atom  potentials  that  have  proven  to  be  successful  in  modeling  a 
large  number  of  organic  crystals  [4, 14].  In  particular,  we  assume  that  (1)  the  intermolecular 
interactions  depend  only  on  the  interatomic  distances;  (2)  the  interaction  potential  can  be 
separated  in  contributions  identified  as  van  der  Waals  and  electrostatic;  and  (3)  the  same  type  of 
van  der  Waals  potential  is  used  for  the  same  type  of  atoms,  independent  of  their  valence  state. 

We  consider  the  case  of  a  crystal  composed  of  rigid  molecules  with  one  molecule  in  the 
asymmetric  unit.  In  this  case,  the  maximum  number  of  degrees  of  freedom  is  12  and  corresponds 
to  the  6  unit  cell  constants  (Ij,  Ij,  I3,  a,  P,  and  y),  three  rotations  (Sj,  62,  and  63),  and  3 
translations  (Tj,  Tj,  and  T3)  of  the  rigid  molecule.  Let  N  be  the  number  of  molecules  per  unit  cell, 
each  molecule  of  index  i  having  nj  atoms,  such  that  a  given  atom  is  identified  by  indices  ia.  The 
position  of  each  unit  cell  in  the  bulk  is  given  by  n=nili+n2l2+n3l3,  where  nj,  n2,  and  n3  are  any 
combination  of  integers  and  Ij  (i- 1,2,3)  are  the  vectors  defining  the  edges  of  the  unit  cell.  The 
reciprocal  space  vectors  k„  are  defined  using  the  matrix  of  basis  vectors  h=[lj  ,12,13] 
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(4) 


K  = 


/  ^ 
nil 


with  nil,  ^3  integers. 

In  the  present  treatment  we  approximate  the  intermolecular  lattice  energy  of  the  crystal  as  a 
pairwise  sum  of  Buckingham  (repulsion  and  dispersion)  and  Coulomb  potentials  of  the  form 


with 


^  ■  ^rep  ^disp"^  ^Coul’ 


^rcp  “  ZEE  AjojjpCXpC  tniajp)’ 


2  n  i=l  «=1  j=l  P=1 


(5) 


(6) 


1  «  N  ®i  N  “j  , 

^disp  “  tE^E  E  E  E  ^ajP^^niajP’ 

2  n  i=l  a=l  j=l  P=1 


(7) 


and 


-  CO  N  N 

^Coul  “  tE*E  E  E  E  Qio^jP'^t'niojjpj 

/  n  i=l  a=l  j=l  p=l 


(8) 


where  r„i„jp=|ri„-rjp+n|  is  the  distance  between  atom  a  of  molecule  i  and  atom  P  of  molecule  j, 
located  in  the  cell  determined  by  the  set  (ni,n2,n3).  The  summation  in  equations  (6)-(8) 

n 

excludes  the  terms  with  H,  when  |  nj  =0. 

In  performing  these  lattice  sums,  special  attention  should  be  paid  to  the  dispersion  and 
Coulombic  sums  (equations  [7]  and  [8],  respectively)  due  to  their  slow  convergence  over  the 
infinite  periodic  lattice.  For  the  treatment  of  these  terms,  we  adopt  the  method  of  accelerated 
convergence  of  die  crystal-lattice  potential  sums  introduced  by  Williams  [15, 16].  In  this 
procedure,  a  general  lattice  sum  Sjj,  of  the  form 
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j  »  N  N 


=  in.  m 

2  n  i=l  a=l  j=l  p-1 


(9) 


is  transformed  by  multiplying  the  terms  of  the  series  by  a  convergence  function  0„(r): 

Z  Z  Z  WjP 

2  n  i=l  a=l  j=l  p=l 


"  |e*z  e  i  e 

2  n  i=l  a=l  j=l,  p=l 


(10) 


The  function  ^ni(r)  rapidly  decreases  to  zero  with  increasing  r  and  ^>in(0)=l. 

As  a  result  of  this  transformation,  the  first  term  in  equation  (10)  converges  much  faster  than 
the  original  sum  S^,  in  equation  (9)  and  can  be  evaluated  in  the  real  space.  The  second  term  in 
equation  (10)  is  still  a  slowly  varying  function  of  r,  but  can  be  determined  at  a  much  higher 
convergence  rate  when  it  is  Fourier  transformed  and  evaluated  in  the  reciprocal  lattice  vector 
space.  The  convergence  function  ^►^(r)  is  chosen  to  be  the  normalized  incomplete  gamma 
function 


^.(r)  = 


— —  J 
r(m/2)  ,4^ 


tm/2-lg-tjjt 


(11) 


of  index  m=l  for  the  electrostatic  case  and  m=6  for  the  dispersion  sums.  The  adjustable 
parameter  t|  in  equation  (11)  has  the  dimension  of  inverse  length  and  determines  the  relative 
contributions  of  the  real-  and  reciprocal-space  terms.  The  contribution  of  reciprocal-space  series 
increases  as  q  increases.  Details  of  the  numerical  manipulations  can  be  found  in  Williams  [15, 
16],  Nijboer  and  Dewette  [17],  and  Karasawa  and  Goddard  [18]. 
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Under  the  assumption  of  neutrality  for  the  total  electrostatic  charge  of  each  molecule,  the 
expression  of  the  Coulombic  energy  determined  using  the  accelerated  convergence  method  can 
be  written  as  superposition  of  a  summation  evaluated  in  the  real  space,  one  evaluated  in  the 
reciprocal  space,  E^,  plus  a  set  of  correction  terms,  E^^,,  that  correspond  to  exclusion  of  the 
self-energy  and  subtraction  of  intramolecular  interactions; 

Ecoui  =  Ei,  +  Ei  +  (12) 


The  analytic  expressions  of  these  sums  are: 


-J  »  N  N  Iij  q  q 

Ei  =  E  Z  E  ^erfc(a,), 

2  n  i=l  a=l  j=l  P=1  Tnijtjp 


(13) 


1 


27tV 


Z  |Si(kJ| 


2  exp(-bi) 


(14) 


and 


E‘  = 

■^cor 


N  “i 

E  I  4 

\i=l  a=l 


^  -I  N  q.  q 

-  ^  I  E  E  ^  erf(a,). 

2  i=l  a=l  p=a+l  r^jjjp 


(15) 


hi  equations  (13)  and  (15),  erfc(...)  and  erf(...)  are  the  complementary  error  and  error  functions, 
respectively.  In  equation  (14)  V=li(l2xl3)  is  the  volume  of  the  unit  cell.  The  parameters  in 
equations  (13)  and  (14)  are  defined  as  ai=riirni„jp  and  bi=iik^T|i.  The  term  Si(kn,)  denotes  the 
structure  factor  defined  as 


N  “i 


Si(kJ  =  E  E  qiaexp(-2Kik^rj^). 

i=l  a=l 


(16) 
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Similarly,  application  of  the  accelerated  convergence  technique  to  the  case  of  the  dispersion 
energy  yields  the  following  summations  [16]: 


and 


^dir 


If  f  2  (1.  .  a//2)exp(-a,^), 


n  i=l  a=l  j=l  p=I  1*° 


niajp 


(17) 


,6  TC 


,9/2 


2h 


•^V  k  #0  \  0^ 


exp(-bg), 


(18) 


6V 


‘  ^  \a 

I  I  ^ 

Vi==l  a=l 


12 


N  “i 

EEC.. 

.  Lj  i-i  '^laio 

Vi=l  a=l 


-4  I  "f  E  [2-(2  •  2a/  *  a/)exp(-a/)]  Ss, 

2  i=l  a=l  p=a+l  rj„jp 


(19) 


The  &st  term  in  equation  (17)  represents  the  contribution  in  the  limit  of  the  Fourier 

transform  of  the  second  term  in  equation  (10).  The  second  term  in  equation  (19)  is  the  self¬ 
energy,  and  the  last  term  of  equation  (19)  represents  the  contribution  of  intermolecular  dispersion 
terms.  The  definitions  of  parameters  ag  and  bg  are  similar  to  those  previously  given  for  Uj  and  bj, 
but  correspond  to  a  different  adjustable  parameter,  qg.  Similarly,  the  definition  of  Sg(k^)  is  that 
of  equation  (16),  but  with  charges  qi„  replaced  by  potential  terms  (Ci„)’'^. 

For  rigid  molecules  (which  are  assumed  in  this  study),  the  terms  corresponding  to  self¬ 
energies  and  intramolecular  energies  do  not  contribute  to  the  derivatives  of  the  lattice  energy.  In 
the  case  when  the  molecules  in  the  unit  cell  are  related  by  symmetry,  the  general  expressions 
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equations  (13)-(15)  and  equations  (17H19)  representing  the  energy  of  a  unit  cell  can  be  further 
reduced.  Full  details  of  this  particular  case  are  given  in  Williams  [16]. 

In  practice,  the  interactions  in  real  space  are  calculated  for  atom  pairs  separated  by  no  more 
than  a  given  cutoff  distance  (r^J.  Interactions  between  atom  pairs  separated  by  distances  larger 
than  r^t  are  assumed  to  be  zero.  Once  r^^,  is  specified,  the  parameters  q,  and  qg  are  then  chosen 
to  ensure  the  proper  convergence  of  the  lattice  sums  in  both  real  and  reciprocal  space  while 
limiting  the  computational  expense  of  the  calculations  of  the  sums  [18, 19]. 

An  optimal  set  of  potential  parameters  can  be  obtained  by  minimizing  a  general  function  of 
the  form 


R  =  E  Wy 

iO 


'  aE' 

ispj 

[apj 

+  w'(E-Ef 


(20) 


at  the  observed  equilibrium  structure.  The  function  R  in  equation  (20)  represents  a  weighted 
superposition  of  forces  and  torques  plus  the  square  of  the  difference  between  the  calculated 
lattice  energy  E  and  the  measured  value  E®  [20].  The  weight  elements  Wy  are  calculated  as 
[h*  VH]:*  ,  where  [H]is  the  Hessian  matrix  (Hij=a^E/apiapj)  and  V  is  a  diagonal  matrix  with 
elements  Vii-o^(pi),  equal  to  the  allowed  error  threshold  of  different  structural  parameters.  The 
static  lattice  energy  E°  in  equation  (20)  can  be  estimated  from  the  experimental  enthalpy  of 
sublimation  by  using  the  relation  [21]:  =  E+J^+2RT,  where  E  is  the  lattice  energy  and 

Kq  is  the  zero-point  energy.  We  have  made  the  approximation  that  the  static  lattice  energy  is 
equal  to  the  measured  sublimation  enthalpy  of  the  RDX  crystal  (- 130.1  kJ/mole)  [22]. 
Discussions  of  alternative  techniques  used  for  optimization  of  potential  parameters  can  be  found 
in  Hsu  and  Williams  [20]. 

Several  of  the  potential  parameters,  including  all  of  the  electrostatic  charges,  were  assigned 
before  optimization  of  equation  (20)  was  attempted.  The  assignment  of  the  electrostatic  charges 
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poses  a  problem  in  that  the  atom-centered  monopole  charge  is  not  an  observable  quantity  and 
cannot  be  obtained  directly  from  either  experiment  or  first  principles  calculations.  There  are 
several  schemes  for  evaluation  of  charges  by  empirical  partition  or  by  using  a  quantum 
mechanically  derived  wave  function  [23].  We  have  chosen  to  make  the  assignments  for  the 
atom-centered  monopole  charges  by  using  the  set  that  best  reproduces  the  quantum  mechanically 
derived  electrostatic  potential  that  is  calculated  over  grid  points  surrounding  the  van  der  Waals 
surface  of  the  RDX  molecule.  This  method  of  fitting  to  the  electrostatic  potential  was  proposed 
by  Breneman  and  Wiberg  [24],  and  is  incorporated  in  the  Gaussian  94  package  of  programs  [25] 
under  the  keyword  CHEIJPG.  This  method  presents  the  advantages  of  having  a  higher  density  of 
points  and  a  better  selection  procedure,  which  ensures  a  decrease  of  rotational  variables  observed 
with  similar  methods.  The  electronic  structure  calculations  were  performed  at  the  MP2/6-31G** 
level  [26-28].  The  atomic  arrangement  of  the  molecule  used  in  these  calculations  was  consistent 
with  the  crystallographic  configuration  [2].  A  total  number  of  1 1,900  points  were  used  in  the  fit, 
leading  to  a  root-mean-square  deviation  of  0.00259.  The  resulting  electrostatic  charges  are  given 
in  Table  1. 


The  number  of  parameters  to  be  fit  was  further  reduced  by  using  traditional  combination  rules 
to  obtain  the  heteroatom  parameters  from  homoatom  parameters: 

■^ap  ^  V'^aa'^pp  » 


®ap  “  (^cMi'^^Pp)/^’ 


--ap  V  ««  PP " 


We  have  used  the  published  6-exp  parameters  [29]  for  the  H-H  and  C-C  nonbonded 
interactions  and  have  optimized  the  values  of  N-N  and  0-0  interactions.  It  has  been  shown  that 
there  is  a  strong  correlation  between  the  nonbonded  potential  parameters  A  and  B  [28]. 


Thus,  we  have  used  published  values  [27]  for  B  parameters  for  the  N-N  and  0-0  interactions 
and  optimized  the  corresponding  A  and  C  parameters. 
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Table  1.  Electrostatic  Charges  for  an  RDX  Molecule 


Atom® 

Charge/  e 

Cl 

-0.092940 

Q 

-0.015018 

Q 

-0.045006 

Ni 

0.038600 

N2 

-0.331673 

N3 

-0.340479 

N4 

0.557379 

N5 

0.781316 

Ne 

0.772447 

Oi 

-0.356222 

O2 

-0.366350 

O3 

-0.369139 

O4 

-0.380884 

O5 

-0.358609 

06 

-0.378756 

Hi 

0.130328 

H2 

0.174188 

H3 

0.161566 

H4 

0.146993 

H5 

0.113532 

E. 

0.158728 

*  The  atom  indices  are  consistent  with  the  labels  in  Figure  2. 
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After  making  these  assignments,  a  two-step  calculation  was  performed  iteratively  until  a 
suitable  parameter  set  was  obtained.  The  first  step  in  the  iteration  involved  selection  of  a  trial 
parameter  set  through  a  combination  of  trial-and-error  adjustments  of  parameters  and 
minimization  of  equation  (20).  The  model  using  the  trial  parameter  set  was  then  subjected  to 
symmetry-constrained  molecular  packing  calculations  using  the  PCK91  set  of  programs  [29] 
(described  in  the  next  section).  These  two  steps  were  repeated  until  a  parameter  set  that 
reproduced  experimental  information  to  the  desired  tolerance  was  obtained.  The  PCK91  set  of 
programs  allows  analytical  derivatives  of  the  potential  only  in  the  direct  space.  To  circumvent 
the  problem  arising  from  the  inability  of  the  PCK91  programs  [29]  to  evaluate  the  derivatives  of 
the  reciprocal  space  terms  (equations  [14]  and  [18]),  we  set  the  cutoff  distance  in  the  real  space 
to  17.0  A.  By  choosing  such  a  large  cutoff,  we  were  able  to  adjust  the  parameters  t|i  and  rig  such 
that  the  reciprocal  sums  are  smaller  than  the  desired  errors  in  the  evaluation  of  equation  (20). 
Thus,  evaluation  of  the  reciprocal  terms  using  these  parameters  became  unnecessary.  The  values 
obtained  for  t|i=0.1861  and  T|6=0.2304  assured  contributions  of  the  reciprocal  sums  smaller  than 
0.01  kJ/mol.  This  large  cutoff  was  also  used  in  the  subsequent  full-minimization-molecular 
packing  and  symmetry-constrained  Monte  Carlo  calculations  described  as  follows. 

The  quality  of  the  fit  was  judged  according  to  two  discrepancy  factors.  The  first  is  a 
structural  shift  factor  of  the  form 


Fj  =  (5OA0)2  +  (10Ax>2  +  X 

j=i 


' 

100— i 
V  y 


■  I(Aa/, 


(22) 


where  A0  is  the  total  rms  rigid-body  rotational  displacement  (in  radians)  after  minimization.  Ax 
the  rms  total  rigid-body  translational  displacement  (in  A),  1;  the  lengths  of  the  edges  of  the  unit 
cell,  and  ttj  the  angles  of  the  unit  cell.  The  second  discrepancy  factor,  F2,  was  calculated  from  Fj 
by  including  the  weighted  difference  between  the  calculated  lattice  energy  and  the  observed 
lattice  energy.  The  threshold  values  are  those  previously  proposed  by  Starr  and  Williams  [31]: 
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unit-cell  edge  of  1%,  cell  angles  and  molecular  rotations  at  0.02  rad,  molecular  translations  of 
0.05  A,  and  heats  of  sublimation  of  0.02  kJ/mol. 

A  (3x3x3)  block  of  unit  cells  of  RDX  with  each  unit  cell  identical  to  the  structure  determined 
by  single-crystal  neutron  diffraction  experiments  [2]  was  used  in  the  parameterization  of  the 
potential  energy  function.  The  refined  structure  of  RDX  at  room  temperature,  1  atm,  belongs  to 
the  orthorhombic  space  group  Pbca  with  Z=8  molecules  per  unit  cell  (see  Figure  1).  The  cell 
parameters  are  li=13.182  A,  12=11.574  A,  and  13=10.709  A,  giving  the  volume  of  the  unit  cell  as 
1633.8557  k\  As  shown  in  Figure  2,  the  molecule  consists  of  three  alternating  N-NOj  and  CHj 
groups  arranged  in  a  six-membered  C-N  puckered  ring. 


Figure  1.  Unit  cell  of  oc-RDX  crystal  using  refined  coordinates  given  in  Choi  and  Prince  121. 
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The  values  of  the  complete  set  of  optimized  intermolecular  parameters  are  given  in  Table  2, 
and  the  values  of  the  crystal  parameters  resulting  from  the  symmetry-constrained  energy 
minimization  using  these  parameters  are  given  in  Table  3.  A  total  translation  of  0.126  A  and  a 
total  rotation  of  1.241  °  away  from  the  experimental  values  for  the  molecule  in  the  asymmetric 
unit  take  place  in  the  minimization  process.  The  corresponding  values  of  the  discrepancy  factors 
are  Fi=0.88  and  F2=0.57. 

3.  DETAILS  OF  THE  CALCULATIONS 

Three  series  of  calculations — ^molecular  packing,  symmetry-constrained  isothermal-isobaric 
Monte  Carlo,  and  isothermal-isobaric  molecular  dynamics  calculations — ^were  performed  to 
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Table  2.  RDX  Atom-Atom  Potential  Parameters 


Pair 

Gop 

Bop 

A-e-p 

(a-|3) 

(kJ/mol  A®) 

(A-1) 

CkJ/mol) 

H-H 

136.3800 

3.74 

9213.510 

C-C 

2439.3459 

3.60 

369726.330 

N-N 

1668.3316 

3.78 

264795.246 

0-0 

1453.3114 

3.96 

290437.820 

assess  the  quality  of  the  potential  energy  function.  Details  of  the  methods  for  each  of  the 
calculations  are  described  herein.  In  each  of  the  calculations,  the  set  of  refined  coordinates  from 
the  neutron-diffraction  study  [2]  was  used  for  the  initial  geometry  of  the  system. 

3,1  Parlfinp.  Molecular  packing  calculations  are  minimizations  of  the  lattice 

energy  with  respect  to  structural  degrees  of  freedom  of  the  crystal.  For  the  particular  case  of 
crystals  with  one  molecule  in  the  asymmetric  unit  occupying  an  arbitrary  position,  the  maximum 
number  of  structural  degrees  of  freedom  is  12.  A  reduced  number  of  structural  degrees  of 
freedom  might  be  involved,  depending  on  the  symmetry  restrictions  of  different  space  groups. 

Assuming  that  the  crystal  energy  is  known  as  a  function  of  the  structural  lattice  parameters, 
the  equilibrium  crystal  configuration  is  determined  by  the  conditions  of  zero  force  and  torques, 
together  with  the  requirement  that  there  is  a  minimum.  The  search  for  such  a  minimum  can  be 
done  using  a  combination  of  steepest-descent  and  Newton-Raphson  procedures  [15, 32]. 

The  Newton-Raphson  method  is  used  to  minimize  the  lattice  energy  for  configurations  close 
to  the  equilibrium,  when  all  eigenvalues  of  the  Hessian  matrix  are  positive  definite.  In  this  case, 
the  step  for  the  iterative  search  for  a  minimum  energy  from  a  trial  configuration  is  determined  as 

Ap  =  -[H(p)r'G(p),  (23) 
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where  G  is  the  vector  of  first  derivatives  and  His  the  matrix  of  second  derivatives  of  the  lattice 

energy  with  respect  to  the  structural  variables.  This  procedure  is  repeated  until  no  sigmficant 
changes  are  observed  in  the  coordinates,  and  the  first  derivatives  of  the  energy  are  close  to  zero. 


In  those  situations  when  the  Hessian  is  not  positive  definite  for  the  trial  configuration,  the 
steepest-descent  method  is  used.  In  this  procedure,  the  searching  step  is  taken  along  the  direction 
of  the  local  downhill  gradient  [32].  This  method  has  a  low  rate  of  convergence  near  the 
minimum,  but  Can  be  efficiently  used  to  determine  a  configuration  for  which  the  Newton- 
Raphson  method  can  be  used. 

The  molecular  packing  program  PCK91  [29]  was  used  in  the  detemunation  of  a  suitable 
parameter  set.  This  program  was  used  to  find  energy  minima  for  systems  described  by  the  trial 
parameter  sets.  This  program  calculates  crystal  lattice  sums  using  the  accelerated  convergence 
method  in  section  2,  and  the  first  and  second  derivatives  of  the  crystal  lattice  energy  are 
evaluated  analytically.  The  space  group  symmetry  is  maintained  throughout  the  energy 
minimization.  This  reduces  the  number  of  independent  variables  in  the  minimization  procedure, 
resulting  in  a  significant  decrease  of  computational  time  when  compared  to  a  nonconstrained 
energy  minimization.  The  crystallographic  parameters  varied  in  the  minimization  using  PCK91 
[29]  are  the  three  dimensions  of  the  unit  cell  and  the  three  respective  translations  and  rotations  of 
a  central  molecule  (from  which  the  positions  of  all  other  molecules  in  the  unit  ceU  can  be 
determined  through  symmetry  operations).  The  angles  of  the  unit  cell  were  firozen  at  90°. 

A  necessary  condition  in  assessing  the  accuracy  of  the  empirical  potential  energy  is  that  the 
structure  corresponding  to  the  energy  minimum  of  the  model  should  maintain  the  observed  space 
group  symmetry  [34].  In  order  to  test  our  potential  for  such  a  requirement,  we  have  used  an 
algorithm  recently  proposed  by  Gibson  and  Scherega  [35]  for  efficient  minirmzation  of  the 
energy  of  a  fully  variable  lattice  composed  by  rigid  molecules.  This  algorithm  makes  use  of 
Gay’s  [36]  secant-type  unconstrained  minimization  solver  (SUMSL)  routine  to  miniimze  the 
lattice  energy.  The  gradients  of  the  energy  with  respect  to  generalized  coordinates  (i.e.,  6Z  rigid- 
body  parameters  of  the  Z  molecules  in  the  unit  cell  and  the  six  lattice  parameters)  are  evaluated 
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analytically.  The  nonbonded  interactions  were  cut  off  at  a  distance  Qo  with  a  cubic  feather 
(spline)  applied  over  the  distance  Po-Qo  to  ensure  that  the  energy  and  its  first  derivatives  are 
continuous  everywhere.  Here,  the  P  and  Q  parameters  specify,  respectively,  the  start  and  the  end 
of  the  feather,  and  a  is  the  value  of  interatomic  potential  at  which  E  =  0  and  dE/dr  <  0.  Long- 
range  Coulombic  interactions  are  evaluated  based  on  an  Ewald  type  of  transformation. 

At  the  begiiming  and  at  the  end  of  lattice  energy  minimization,  the  symmetry  operations  that 
transform  the  molecules  in  the  unit  cell  are  computed,  and  if  these  operations  remain  unchanged, 
the  space  group  is  considered  conserved;  otherwise,  a  new  space  group  is  deduced  based  on  the 
final  symmetry  relations.  Some  of  the  symmetry  transformations  can  be  lost  during  energy 
minimization  but  are  regained  before  convergence  takes  place  [34]. 

We  have  used  the  algorithm  as  implemented  in  the  program  LMIN  [37]  to  analyze  the 
accuracy  of  the  proposed  intermolecular  potential  for  the  RDX  crystal.  In  all  of  these 
calculations,  the  parameter  P  was  set  equal  to  ~Q-5. 


3.2  Symmetry-Constrained  NPT  Monte  Carlo  Calculations.  The  intermolecular  potential 
was  further  evaluated  based  on  its  performance  in  Monte  Carlo  calculations  in  the  NPT  ensemble 
[19].  These  calculations  were  used  to  determine  some  of  the  crystal  parameters  as  functions  of 
temperature  and  pressure  as  predicted  by  the  model.  In  these  calculations,  we  have  constrained 
the  crystal  symmetry  as  in  the  first  series  of  molecular  packing  calculations  using  PCK91  [29]. 
The  simulation  box  contains  all  unit  cells  from  -3,3  along  the  Ij  axis,  -3,4  along  the  I2  axis,  and 
-3,4  along  the  I3  axis,  with  the  central  cell  occupying  position  (0,0,0).  The  positions  of  the 
molecules  inside  the  unit  cell,  considered  as  rigid  entities,  are  determined  from  the  coordinates  of 
a  single  molecule  in  the  unit  cell,  which  is  denoted  hereafter  as  (T).  The  coordinates  of  all 
remaining  molecules  in  the  crystal  of  RDX  are  obtained  from  the  fractional  positions  of  (T)  by 
applying  the  operators  of  the  crystal  space  group  Pbca  [38].  The  position  of  (T)  inside  the  unit 
cell  is  described  by  a  set  of  fractional  coordinates  of  the  mass  center,  and  the  molecular 
orientation  is  described  by  a  set  of  quaternions  [19, 39].  Because  of  the  symmetry  constraints 
within  the  simulation,  each  molecule  in  this  system  has  the  same  number  and  spatial  arrangement 


of  neighbors  as  (T)  at  each  step  in  the  NPT-MC  simulation.  Using  these  symmetry  constraints, 
the  total  potential  energy  of  a  system. 


,  «  N  ”i  N  “j 

E  •  i  E*  E  I  I  I 

2  n  i=l  a=l  j=l  P=1 


(24) 


can  be  written  in  terms  of  the  potential  energy  of  a  single  molecule  (T), 


XI  “  “t  N  “j 

E  =  ^  E*  E  E  E 

2  n  a=l  j=l  |J=1 


(25) 


Therefore,  in  this  type  of  simulation,  it  is  necessary  to  calculate  only  the  potential  energy  for  the 
single  molecule  (T)  (equation  [25])  rather  than  perform  the  summations  over  all  molecules  as 
given  in  equation  (24).  The  computational  demands  of  symmetry-constrained  NPT-MC 
simulations  are  much  less  than  nonconstrained  NPT-MC  simulations,  since  the  phase  space  that 
is  sampled  is  reduced  to  the  structural  parameters  corresponding  to  a  single  molecule.  We  found 
in  this  and  a  previous  study  [40]  that  the  sunulations  provide  averages  that  do  not  differ  radically 
from  non-s)rmmetry-constrained  NPT-MC  predictions.  This  allows  for  rapid  testing  of 
intermolecular  interaction  potentials  using  the  method  of  Monte  Carlo.  Poor  potentials  can  be 
eliminated  quickly  with  these  inexpensive  calculations,  while  potentials  that  show  reasonably 
good  results  in  these  simulations  can  be  subjected  to  more  rigorous  nonconstrained  simulations 
such  as  those  described  later. 

The  general  method  followed  in  performing  the  MC-NPT  calculations  is  that  described  in 
Allen  and  Tildesley  [19].  Here  we  will  discuss  only  the  main  steps  of  this  procedure.  The 
sampling  in  the  MC  walk  was  performed  over  the  three  cell  constants  corresponding  to  unit  cell 
edges,  the  fractional  coordinates  of  the  mass-center,  and  the  set  of  quaternions  of  (T).  The  angles 
of  the  unit  cell  were  frozen  at  90°.  At  every  step  during  the  Markov  walk,  a  trial  to  randomly 
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modify  both  the  position  and  orientation  of  (T)  and  one  edge  of  the  unit  cell  was  attempted. 
These  types  of  displacements  were  performed  according  to  the  relations: 


sr  =  i=x,y,z 

qr  =  i=l,2,3.4 

l_new  i.12,3,  (26) 

where  and  ^3  are  random  numbers  uniformly  distributed  in  the  interval  (0,1),  and  As^, 
Aq^„,  and  Al„„  are  the  maximum  allowed  changes  of  the  fractional  coordinates  of  the  mass 
center,  quaternions,  and  unit  cell  dimensions,  respectively.  The  maximum  displacements  are 
continuously  adjusted  using  the  procedure  described  in  Allen  and  Tildesley  [19]  to  maintain  an 
acceptance/rejection  ratio  of  about  50%. 

After  each  trial  move,  the  new  coordinates  of  the  molecules  in  the  crystal,  the  new  volume  V, 
and  the  new  interaction  energy  (equation  [25])  were  calculated.  The  acceptance  or  rejection 
of  each  trial  move  was  decided  according  to  the  probability  min(l,exp[-pW]),  where 

W  =  -Pl(E„.  -  E„„)  +  P(V^  -V„„)]  ♦  N  (27) 

and  P=l/(kT).  The  move  was  accepted  if  W  <  0.  If  W  >  0,  a  new  random  number  C  distributed 
in  the  interval  (0,1)  was  generated  and  compared  to  exp(-pW).  The  move  was  accepted  if 
C  <  exp(-pw)  and  rejected  otherwise.  Averages  of  different  quantities  such  as  the  lattice  energy, 
lattice  dimensions,  density,  and  the  positions  of  the  (T)  are  then  determined.  The  fluctuations  of 
a  property  A,  defined  as 


F  =  )/  A'  -  (A)2, 


(28) 


will  be  reported  with  the  averages  for  these  simulations,  as  well  as  for  the  NPT-MD  simulations 
described  in  succeeding  text.  In  all  the  MC  simulations,  warm-up  walks  of  1,000-3,000  MC 
steps  were  done.  The  Markov  walk  was  performed  for  a  total  number  of  steps  varying  from 
100,000  at  T=4.2  K  to  250,000  at  300  K. 

A  problem  of  interest  in  NPT  simulations  is  the  evaluation  of  the  average  pressure.  For  a 
system  of  rigid  polyatomic  molecules.  Smith  [41 ,  42]  has  shown  that  the  pressure  can  be 
calculated  using  the  molecular  virial  theorem  if  the  Hamiltonian  of  the  system  is  rewritten  in 
terms  of  the  scaled  coordinates  of  the  molecular  center  of  mass  (s—r/V^^)  and  of  the 
intramolecular  distances  relative  to  the  center  of  mass  (di„=rj„-rj),  invariant  at  volume  scaling. 

A  more  general  approach  has  been  introduced  by  Parrinello  and  Rahman  [43],  who  have  used  the 
full  microscopic  stress  tensor  for  an  N-particle  system  in  a  periodically  repeating  MD  cell.  Nose 
and  Klein  [44]  have  extended  the  calculation  of  the  stress  tensor  to  include  the  long-range 
Coulombic  interactions  and  have  applied  this  to  rigid  polyatomic  molecules.  We  present  in  the 
Appendix  a  generalization  of  the  treatment  based  on  the  accelerated  convergence  technique  given 
by  Karasawa  and  Goddard  [18],  to  formulate  the  stress  tensor  for  both  the  Coulombic  and 
dispersion  potentials  assuming  a  system  of  rigid  molecules  in  a  periodic  lattice.  The  calculated 
pressure,  as  given  in  equation  (A9),  has  been  used  to  provide  an  independent  verification  of  the 
reliability  of  the  codes  developed  in  this  study. 

3.3  NPT  Molecular  nvnamics  Calculations.  The  intermolecular  interaction  potential  for 
RDX  was  subjected  to  a  final  and  more  rigorous  test:  NPT  molecular  dynamics  simulations,  in 
which  there  are  no  constraints  other  than  the  assumption  of  rigid-body  molecules.  Single 
trajectories  at  4.2, 100, 200, 250, 273.15, 300,  and  325  K  and  zero  pressure  were  integrated  for  at 
least  10  ps.  The  calculations  were  performed  using  the  MDSCPC4  set  of  programs  [45].  In  this 
procedure,  the  equations  of  motion  for  both  the  molecules  and  the  simulation  cell  are  integrated 
using  a  fifth-order  Gear  predictor-corrector  [46].  The  rotational  motion  is  handled  using  the 
Evans  quaternion  algorithm  [47]  with  a  fifth-order  Gear  integrator.  The  algorithms  used  are 
described  in  detail  in  Nose^  and  Klein  [44].  A  fixed  time  step  of  2.0x10  s  was  used.  The 
system  consists  of  a  box  containing  27  unit  cells  [a  3x3x3  box  of  unit  cells].  At  the  beginning  of 
each  simulation,  the  unit  cells  are  identical  to  the  experimental  structure  at  T=300  K,  1  atm  [2]. 
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The  initial  velocities  of  the  centers-of-mass  of  the  molecules  were  selected  at  random,  but  subject 
to  the  following  conditions:  the  translational  and  rotational  velocities  were  selected  such  that  the 
simulation  cell  had  no  bulk  motion.  Also  the  translational  and  rotational  velocities  were  scaled 
to  yield  the  desired  temperature.  The  system  was  first  integrated  for  2,000  steps  (4  ps)  under 
conditions  of  constant  volume  in  order  to  relax  the  system  away  from  the  initial  state.  During 
this  period,  the  velocities  were  scaled  at  every  five  steps  in  order  that  the  internal  temperature  of 
crystal  mimic  the  imposed  temperature.  And  additional  3,000  steps  (6  ps)  were  then  integrated 
(at  constant  volume)  without  additional  scaling  of  the  velocities.  Finally,  the  isothermal-isobaric 
trajectory  was  integrated  for  at  least  5,000  steps  (10  ps)  for  each  temperature  at  zero  pressure. 
During  the  first  2,000  steps  (4  ps)  of  each  NPT  trajectory,  the  velocities  were  scaled  at  every  5 
steps  as  described  previously.  Averages  were  obtained  from  properties  calculated  at  subsequent 
integration  steps  in  the  NPT  simulation. 

Summations  given  in  equations  (13)-(19)  were  adapted  for  inclusion  of  minimum-image 
periodic  boundary  conditions  in  all  dimensions  in  these  simulations  [19].  The  interactions  are 
determined  between  the  sites  (atoms)  in  the  simulation  box  and  the  nearest-image  sites  within  the 
cutoff  distance.  In  these  calculations,  the  sum  over  reciprocal  vectors  was  included  for  the 
Coulombic  interactions,  subject  to  the  summation  limit  given  by  (mj^  +  m2^-H  m3^)  ^12^  (see 
equation  [4]).  Long-range  corrections  for  the  potential  energy  and  virial  contribution  due  to 
dispersion  were  calculated  using  the  standard  techniques  [19].  The  cutoff  distance  here  is  80% 
of  one-half  the  shortest  perpendicular  distance  between  two  faces  of  the  simulation  cell, 
assuming  unit  cell  sizes  measured  at  300  K,  1  atm  [2].  This  corresponds  to  a  cutoff  distance  of 
12.8508  A.  We  monitored  the  simulation  box  size  to  be  certain  that  the  shortest  perpendicular 
distance  between  two  opposite  faces  of  the  box  never  became  less  than  twice  the  cutoff  distance 
to  avoid  violating  the  minimum  image  condition  during  the  simulation  [19]. 

4.  RESULTS  AND  DISCUSSIONS 

4.1  Molecular  Packing  Calculations.  The  results  of  the  molecular  packing  calculations  are 
given  in  Table  3.  In  the  full  minimization  of  the  RDX  crystal  energy,  the  space  symmetry  is 
maintained  with  less  than  0.003°  deviations  of  the  lattice  angles  from  experiment.  Also,  the 
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maximum  deviation  of  the  fractional  coordinates  of  the  centroids  [35]  is  0.006.  The  deviation  of 
the  Euler  angles  from  experiment  is  no  greater  than  1.183°. 

Since  accelerated  convergence  is  not  used  for  evaluation  of  the  1/r®  attractive  terms  in  the 
nonbonded  potential  in  the  full  minimization  of  the  RDX  crystal  energy,  dependence  of  the 
lattice  energy  on  the  cutoff  distance  is  expected.  This  is  shown  (see  Table  3)  by  results  obtained 
for  different  Q  parameters.  For  values  of  Q  larger  than  10,  a  small  influence  on  the  total  energy 
per  molecule  and  almost  no  effect  on  geometrical  parameters  of  the  lattice  are  found.  In  addition, 
these  results  are  in  very  good  agreement  with  those  obtained  by  using  the  symmetry  constraints 
procedure  and  accelerated  convergence  technique. 

4.2  Svmmfttrv-rnnstrained  NPT  Monte  Carlo  Calculations.  Ideally,  the  structure  predicted 
in  a  crystal-packing  calculation  should  reproduce  the  crystal  structure  at  temperatures  close  to  0 
K.  Although  experimental  data  in  this  temperature  region  or  extrapolated  values  from 
experiments  performed  at  much  higher  temperatures  are  not  available,  we  have  performed  NPT- 
MC  calculations  at  4.2  K  to  provide  this  comparison.  Additionally,  we  investigated  the  influence 
of  thermal  effects  on  the  crystallographic  stmcture  by  performing  MC  calculations  in  the  NPT 
ensemble  at  temperatures  of  T=100, 200,  and  300  K  and  pressures  of  1  atm  and  500  atm.  The 
averages  and  fluctuations  corresponding  to  lattice  constants,  lattice  potential  energy,  density,  and 
the  positions  and  orientation  of  the  molecule  in  the  asymmetric  unit  are  given  in  Table  4. 

As  evident  in  Table  4,  the  effect  caused  by  the  increase  in  temperature  from  4.2  K  to  300  K 
at  1  atm  is  the  expansion  of  the  lattice  dimensions.  The  most  significant  increase  of  0.53  A  takes 
place  in  a  direction  parallel  to  the  I3  axis  over  this  temperature  range.  Correspondingly,  the 
overall  density  of  the  crystal  decreases  from  1.796  g/cm^  at  4.2  K  to  1.689  g/cm^  at  T=300  K.  A 
similar  trend  is  observed  for  calculations  at  P=500  atm.  The  unit  cell  dimensions  are  within  4% 
of  the  experimental  values  at  T=300  K  and  1  atm,  with  the  largest  deviation  being  the  I3 
parameter,  which  is  0.43  A  larger  than  the  experimental  value. 
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.  NPT  Monte  Carlo  Averages 


No  significant  deviations  of  the  position  and  orientation  of  molecule  (T)  from  the  experimental 
values  were  observed  for  this  range  of  temperatures  and  pressures.  We  found  that  the  maximum 
deviations  of  the  average  values  of  the  fractional  coordinates  and  Euler  angles  from  the  experimental 
data  are  smaller  than  0.006  and  2°,  respectively. 

4,3  NPT  Molecular  Dynamics  Calculations.  Averages  and  fluctuations  of  crystal  structure 
information  calculated  from  the  trajectories  integrated  at  T=4.2  and  300  K  are  given  in  Table  5. 
These  averages  include  lattice  dimensions,  mass-center  fractionals,  and  Euler  angles  of  the  eight 
molecules  of  the  unit  cell.  The  averages  for  each  of  the  eight  molecules  in  the  unit  cell  are  over  the 
27  unit  cells  within  the  simulation  box.  The  results  at  4.2  K  are  in  close  agreement  with  the 
molecular  packing  and  symmetry-constrained  NPT-MC  results. 

The  time  histories  of  the  lattice  parameters  (Ij,  l2»  Is?  P,  and  y)  and  the  volume,  pressure,  and 
rotational  and  translational  temperatures  for  the  T=300  K  trajectory  are  given  in  Figures  3  and  4. 
The  trajectory  is  well  behaved;  each  property  oscillates  about  the  average  value  given  in  Table  5  for 
the  duration  of  the  trajectory.  The  lattice  dimensions  deviate  no  more  than  2%  from  the 
experimental  values.  The  I3  prediction  is  in  much  better  agreement  with  experiment  (within  0.02  A) 
than  the  predictions  of  the  NPT-MC  calculations.  Figure  5  provides  a  comparison  with  experiment 
[2]  of  the  averages  of  orientational  parameters  of  the  eight  molecules  in  the  unit  cell.  It  is  evident 
that  rotational  disorder  is  small  and  the  center-of-mass  positions  of  the  molecules  are  very  close  to 
those  determined  experimentally  [2]. 

Trajectories  at  T=100, 200, 250, 273.15, 293,  and  325  K  were  also  integrated  for  10  ps  to  provide 
information  on  the  behavior  of  the  system  with  increasing  temperature.  Averages  of  the  lattice 
dimensions  resulting  from  these  simulations  are  given  in  Table  6.  This  information  can  be  used  to 
extract  the  thermal  expansion  coefficients  by  using  the  relations: 
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Fluctuations  [Eq.  (28)]  of  MD-NPT  predictions  in  parentheses. 
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Fluctuations  [Eq.  (28)]  of  MD-NPT  predictions  in  parentheses. 
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TIME  (ps) 


Figures.  Time  histories  of  lattice  parameters  (li.  1-,.  L.  a.  (3.  y)  for  isoihermal-isobaric 
trajectory  corresponding  to  T-^300  K.  0  atm. 
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4.  Time  histories  of  unit  cell  volume,  pressure,  rotational  temperature  (TFRotl’). 

and  center-of-mass  translational  temperature  (Trtransli  for  isothermal-isnharic 

trajectory  corresponding  to  T=300  K.  0  atm. 
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nparison  of  time-averaged  center-of-mass  fractional  positions 
Convention)  with  experiment  at  T=300  K.  0  atm. 


Table  6,  NPT-MD  Lattice  Dimensions  vs.  Temperature 


Lattice  Dimensions 


It 

I2 

h 

Volume 

T(K) 

(A) 

U 

(A) 

(A) 

(A') 

Calculated 

Calculated 

Calculated 

Calculated 

4.2 

13.290 

11.654 

10.610 

1643.294 

100 

13.321 

11.699 

10.646 

1659.098 

200 

13.360 

11.748 

10.683 

1676.732 

250 

13.385 

11.771 

10.710 

1687.412 

273.15 

13.389 

11.788 

10.726 

1692.880 

293 

13.394 

11.796 

10.731 

1695.451 

300 

13.396 

11.798 

10.732 

1696.150 

325 

13.416 

11.815 

10.750 

1703.983 

Coefficients  of  Fit  of  Lattice  Parameters  to  Quadratic  in  Temperature 


Parameter^ 

It 

^2 

I3 

Volume 

A 

6.307055  X  10-^ 

-1.512549X  10*’ 

8.828888  x  lO’’ 

9.259430  X  10*^ 

B 

-3.229344  X  lO'^ 

6.484753  x  10"* 

-1.068887  X  10-5 

-3.228835  x  lO'* 

C 

13.79879 

11.61935 

10.65897 

1710.723 

*  A,  B,  and  C  in  units  of  A/T^,  A/T,  and  A,  respectively  for  b,  i-1, 2,  and  3. 
A,  B,  and  C  in  units  of  A^/T^,  A^/T,  and  A^,  respectively  for  the  volume. 


where  a  is  the  coefficient  of  linear  expansion  of  the  material  and  X  denotes  the  length  of  one  of 
the  sides  of  the  unit  cell.  The  coefficient  of  volume  expansion,  P,  has  a  similar  form: 


(30) 
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A  single  value  of  the  linear  expansion  coefficient  at  T=293  K  (63.6  x  10‘®  K  has  been  reported 
[49],  which  implies  that  thermal  expansion  is  the  same  in  each  of  the  three  dimensions.  Two 
volume  expansion  coefficients,  191  x  lO'®  at  293  K  and  250  x  lO'®  K'*  for  the  temperature 
range  293-373  K,  are  also  given.  Figure  6  shows  the  unit-cell  lattice  parameters  and  volume  as 
functions  of  temperature.  The  squares  denote  the  time-averaged  values  obtained  from  the 
trajectories  (listed  in  Table  6),  and  the  lines  denote  fits  of  these  averages  to  quadratic  functions  of 
temperature  over  the  range  250  to  325  K.  The  coefficients  of  each  fit  are  given  in  Table  6.  The 
fits  can  be  used  to  evaluate  equations  (29)  and  (30)  to  provide  the  linear  and  volume  expansion 
coefficients  for  this  model.  The  linear  expansion  coefficients  for  lattice  dimensions  Ij,  I2,  and  I3 
are  34.8  x  10'^  47.5  x  10'®,  and  47.2  x  10‘®  Kr\  which  are  45, 25,  and  26%  smaller  than  the 
experimental  values,  respectively.  The  results  of  the  simulations  differ  from  the  experimental 
values  in  that  the  value  of  the  linear  expansion  coefficient  for  cell  parameter  Ij  is  smaller  than 
those  for  cell  parameters  I2  and  I3,  which  are  approximately  the  same.  The  calculated  value  of 
the  volume  expansion  coefficient  (129.6  x  10  ^  K  *)  at  293  K  is  32%  smaller  than  the  measured 
value  [49]. 

5.  SUMMARY  AND  CONCLUSIONS 

In  this  paper,  we  have  developed  an  intermolecular  potential  for  the  RDX  crystal  based  on 
6-exp  Buckingham  potentials  terms  plus  Coulombic  interactions.  Electrostatic  charges  of 
different  atoms  in  the  RDX  molecule  were  determined  from  fits  to  ab  initio  electrostatic 
potentials  calculated  at  the  MP2/6-31G**  level.  Values  for  heteroatom  potential  parameters 
were  obtained  from  those  for  homoatom  parameters  by  using  traditional  combination  rales. 
Values  for  several  of  the  homoatom  parameters  were  taken  from  published  data  [10, 30], 
including  the  H-H  and  C-C  potential  terms.  The  published  value  of  the  B  parameter  of  the  6-exp 
Buckingham  potential  for  the  N-N  and  0-0  Buckingham  terms  was  used,  and  the  values  of  the 
remaining  A  and  C  parameters  for  these  homonuclear  interactions  were  optimized  by  nonlinear 
least-squares  fitting  of  the  potential  to  a  general  function  that  represents  a  weighted  superposition 
of  forces  and  torques  plus  the  square  of  the  differences  between  the  predicted  and  experimental 
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lattice  energies  for  the  experimentally  determined  crystal  structure.  The  potential  parameters 
were  adjusted  to  reproduce  the  experimentally  determined  space  group  symmetry,  the  lattice 
dimensions  for  the  RDX  crystal,  and  the  enthalpy  of  sublimation. 

Accurate  values  of  the  crystal  lattice  energy  have  been  obtained  in  symmetry-constrained 
molecular  packing  and  NPT-MC  calculations  by  employing  the  technique  of  accelerated 
convergence  for  the  dispersion  and  Coulombic  lattice  sums.  Molecular  packing  calculations  with 
and  without  symmetry  constraints  indicate  very  good  agreement  with  experimental  geometrical 
and  thermochemical  data. 

The  temperature  dependencies  of  the  physical  parameters  of  the  lattice  have  been  investigated 
by  performing  symmetry-constrained  NPT-MC  calculations  in  the  range  4.2—300  K  and  at 
pressures  of  1  and  500  atm.  The  average  values  of  the  physical  lattice  parameters  calculated 
from  Markov  walks  with  100,000-250,000  steps  indicate  that  the  major  modification  of  the 
lattice  takes  place  along  the  I3  axis  with  increasing  temperature.  Molecular  reorientation  and 
translation  of  the  mass-center  fractionals  are  no  greater  than  2°  and  0.006  from  the  experimental 
values,  respectively,  over  the  temperature  and  pressure  ranges  considered.  Lattice  parameters 
calculated  at  T=300  K  and  1  atm  are  within  4%  of  experimental  values. 

A  final  test  of  the  interaction  potential  was  performed  using  isothermal-isobaric  molecular 
dynamics  simulations  at  0  pressure  over  the  temperature  range  4.2—325  K.  The  results  of  these 
calculations  indicate  that  this  model  reproduces  to  within  2%  die  measured  cell  dimensions  of 
the  RDX  crystal  at  300  K.  Additionally,  little  rotational  or  translational  disorder  occurred  in  the 
thermal,  unconstrained  trajectories. 

Linear  and  volume  expansion  coefficients  at  293  K  were  calculated  using  time  averages  over 
the  temperature  range  250-325  K.  The  calculations  differ  from  experiment  in  that  one  of  the 
dimensions  of  the  unit  cells  does  not  have  the  same  linear  expansion  coefficient  as  the  other  two 
dimensions,  as  implied  by  the  single  reported  value  [49]  of  the  linear  expansion  coefficient.  The 
linear  expansion  coefficients  for  the  three  dimensions  are  45, 25,  and  26%  smaller  than  the 
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experimental  value.  The  volume  expansion  coefficient  is  32%  smaller  than  the  measured  value 
[49]. 

These  results  demonstrate  that  the  proposed  potential  provides  an  accurate  description  of  the 
a-form  of  the  RDX  crystal.  The  structure  resulting  from  NPT-MD  simulations  at  300  K  and  0 
pressure  is  in  excellent  agreement  with  measured  structural  information.  The  lack  of  rotational 
and  translational  disorder  of  the  molecules  seen  in  all  of  the  simulations  indicates  that  the 
interaction  potential  has  the  correct  anisotropies  of  the  intermolecular  interactions.  The  good 
agreement  of  the  thermal  expansion  coefficients  with  experiment  at  293  K  suggests  that  the 
model  behaves  properly  over  a  large  temperature  range,  even  though  thermal  information  was  not 
used  in  the  determination  of  the  potential  energy  parameters. 

This  model  appears  to  be  useful  for  prediction  of  nonreactive  processes  occurring  in  an 
a-RDX  crystal.  Refinement  of  this  model  can  be  made  by  including  the  effects  of  intramolecular 
motions,  particularly  of  low-frequency  torsional  motions  of  the  nitro  groups  and  the  ring.  In 
addition,  the  model  can  be  extended  to  reproduce  not  only  geometrical  and  energetic  parameters, 
but  also  spectroscopic  data  for  the  RDX  lattice.  Extensions  of  the  model  to  include  the 
intramolecular  degrees  of  freedom  are  expected  to  facilitate  full  atomistic  investigations  of  the 
dynamics  of  this  energetic  material.  Additionally,  preliminary  results  indicate  that  the 
parameters  determined  for  this  potential  energy  function  are  transferable  to  other  important  cyclic 
nitramines. 
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APPENDIX: 

FORMULATION  OF  THE  STRESS  TENSOR 
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It  has  been  previously  shown  that  for  a  system  of  rigid  molecules  in  a  periodic  system  under 
external  pressure  that  the  internal  stress  tensor  11  can  be  defined  as* 


=  — 
V 


N 


I  (hs,),  + 

i=l  ’I 


(Al) 


where  O  is  the  potential  energy  between  the  molecules  in  the  central  cell  and  the  molecules  in  the 
rest  of  the  crystal.  Here  we  have  assumed  the  Cartesian  coordinates  of  the  atom  a  belonging  to 
molecule  of  index  i  are  represented  as  functions  of  the  molecular  mass  center,  hs  j ,  and  the 

intramolecular  distance,  dj„: 


‘ia 


hsj+d.„, 


(A2) 


where  and  Sj  are  the  firactional  vectors  of  the  atom  (i,a)  and  of  the  molecule  i,  respectively,  and 
h  =  [Ij,l243]  is  the  matrix  of  basis  vectors. 


In  this  appendix  we  give  the  expressions  of  the  Coulombic  and  dispersion  stress  tensor. 

These  results  are  obtained  following  the  method  originally  introduced  by  Nose'  and  Klein^  for  the 
case  of  Coulombic  interactions.  The  corresponding  contribution  of  the  repulsive  energy 
(equation  [6]  in  the  report)  to  the  stress  tensor  is  calculated  only  in  the  real  space  as 


AjttjpBj|^jpexp(  Bj^jprjjjjjjp)(rjji3(jp)g(rjjjj)^. 

2  n  i=l  a=l  j=l  P=1 


N  “i 
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tepey 


(A3) 


*  Nose',  S.,  and  M.  L.  Klein.  Molecular  Physics.  Vol.  50,  p.  1055, 1983. 
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Using  the  identities  =  Vh^j  and  the  Coulombic  stress  tensor  is 

found  to  be  of  the  form  11^  =  +  II^,  where 


1  »  N  “i  N  Dj  /j.  \  (j.  \ 

vnL.,  =  4z*E  1 1  Z  q,.qjp  “7 

2"  i-i-ij-ip-i 
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erfc(a,)  +  — exp(-a5 
Mr 


(A4) 
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1  y  exp(-bt) 
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Si(kJP 


^ev-2^^ - 


-1*1  I  qia2^i(km)c(dia)y[Si(kJexp(-2Tiik^rj^) 

n  i=l  a=l  '■ 


-Sj(-k^)exp(27tik^rjj] 


(A5) 


For  flexible  molecules,  an  additional  term  representing  the  contribution  of  the  potential  term 
should  be  included  (see,  for  example,  Essman  et  al.*). 


Similarly,  the  stress  terms  determined  by  the  dispersion  potential  are  calculated  as: 


^  Essmann,  U.,  L.  Perera,  M.  L.  Berkowitz,  T.  Darden,  H.  Lee,  and  L.  G.  Pedersen.  Journal  of  Chemical  Physics. 
Vol.  103,  p.  8577. 


44 


vnl 


dire  Y 


E  I  I  C^,ju  +  (A6) 

2.  i-i  p-i 


r9/2 


Vltev  =  ^  E  |S«(kjr 

3V  k  ^0 


V^erfc(bg)  +  (-^-  -  -^)exp(-b6^) 


2bi  b. 


CY 


3k^[v/iterfc(bg)-exp(-b6^)/b](kJ^(kJ^ 


q/2  CO  N  Hj 

3V  n  i=l  a=l 


^edcQo^)  + 


2h!  h 


X  I  (Si(k Jexp(-27rik^rj„)  -Sj  (-kj  exp (2nik^rj„)). 


(A7) 


and 


^  -^coreY 


6V 


N  “i 


L  I(Ci«„) 


1/2 


i=l  a=l 


'ey 


(A8) 


In  this  case,  there  is  a  nonzero  contribution  of  the  potential  term  E^^j,  to  the  stress  tensor,  due  to 
the  volume  dependence  of  the  first  term  in  equation  (19).  In  the  previously  listed  expressions, 
indices  oc  and  P  denote  different  molecules  of  the  ensemble,  while  indices  8  and  y  represent  the 
X,  y,  and  z  components  of  the  tensor,  and  6  is  the  Dirac  delta  function. 

Under  hydrostatic  conditions,  at  equilibrium,  the  pressure  in  the  system  can  be  calculated  as 
1/3  of  the  average  of  the  trace  of  the  total  stress  tensor: 
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(A9) 


P 


where  the  total  contributions  from  repulsive,  dispersion,  and  Coulombic  potential  terms  are 
included  in  the  stress  tensorial  terms 


Using  the  identity  SV/dh^^  =  Vh^"g  with  the  general  definition  (Al),  it  follows  that  at 
equilibrium  and  under  hydrostatic  conditions,  an  alternative  formulation  of  the  pressure  is 


P 


J_  /v 

3V  \h  m,  av 


(AlO) 


This  expression  has  been  previously  used  by  Smith^-  for  calculating  the  pressure  in  a  molecular 
system  that  can  expand  or  contract  isotropically. 


^  Smith,  W.  Information  Quarterly  for  Computer  Simtilation  of  Condensed  Phases:  The  CCP5  Newsletter, 
p.  43, 1987. 


^  Smith,  W.  Infonnation  Quarterly  for  Computer  Simulation  of  Condensed  Phases:  The  CCP5  Newsletter, 
p.  14, 1993. 
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